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Abstract 

We present new results on the statistical hadronization of heavy quarks at SPS, 
RHIC and LHC energies. Several new aspects are considered, among them a sepa- 
ration of the collision geometry into a "core" and a "corona" part and an estimate 
of the annihilation rate of charm quark in a hot plasma, together with a critical 
assessment of its influence on the results. For RHIC energies we investigate the 
centrality dependence of J ftp production focusing on the model results for different 
values of the charm production cross section, including its theoretical and experi- 
mental uncertainty. We also study, within this model, the rapidity dependence of the 
J ftp yield. Recent RHIC data from the PHENIX experiment are well reproduced. 
At LHC energy, we update our model predictions for the centrality dependence 
of the J/ip yield and investigate as well the rapidity dependence. We also discuss 
the transverse momentum distributions of J/ip mesons expected from the model 
and provide predictions for a range of values of the expansion velocity at chemical 
freeze-out. Finally, we extend the model to predict T yields in Pb+Pb collisions at 
LHC energy. 



1 Introduction 



The idea of statistical hadronization of charm quarks in nucleus-nucleus collisions [T] has 
led to a series of investigations of J/ip production based on this approach [2|3lH|5|6] . 
Initial interest focussed on the available SPS data for J/ip production in Pb-Pb collisions, 
but the trends for RHIC and LHC were also investigated. An independent approach, 
based on a kinetic model [7 f 8 | 5] has been developed in parallel. Recently the statistical 
hadronization model was extended to the production of T mesons [10J and multiply heavy- 
flavored hadrons 

We have shown earlier [Tf5] that the J/ij) data at SPS energy can be described within the 
statistical approach, but only when assuming that the charm production cross section is 



enhanced by about a factor of 3 beyond the perturbative QCD (pQCD) predictions. For 
RHIC energy, the predictions of our model were shown [5] to be in good agreement with 
the early PHENIX J/if> data [T3] as well as with the relative yields of open charm mesons 
measured by STAR [T5] . 

We focus in this paper on the production of charmonia in nucleus-nucleus collisions at 
RHIC (Au-Au, v /s^7=200 GeV) and LHC (Pb-Pb, v /^7=5.5 TeV) energies within the 
framework of the statistical hadronization model, SHM [T1I5] . Particular emphasis is placed 
on the rapidity and transverse momentum dependence of J/tp and, for LHC energy, also 
T production. Section 2 contains a brief summary of the main ingredients of the model 
along with a description of recent updates concerning the chemical freeze-out parameters, 
the spectrum of open charm mesons and baryons, and the values of charm production 
cross section. Furthermore, for a more realistic description of the centrality dependence in 
nucleus-nucleus collisions, we separate the collision geometry into a 'core' and a 'corona' 
part. 

In Section 3, we provide an estimate of the magnitude of cc annihilation in the expanding 
quark-gluon plasma (QGP). 

In Section 4, we use the model to predict phase space distributions of quarkonia pro- 
duced in nucleus-nucleus collisions. We compare the results with SPS data and with data 
available at RHIC [H] and provide predictions for the LHC energy, where data [15] are 
expected in about two years. In particular the situation at SPS energy is significantly 
modified when taking into account the realistic reaction geometry. We extend the model 
for the calculation of T yields in Pb+Pb collisions at LHC energy. In this case, despite 
the similarity in production rates of bottom at LHC with charm at RHIC [16J, the ap- 
plicability of our model may be questionable, as the thermalization of the bottom quarks 
may be debatable. We discuss below the main sources of uncertainty in these calculations. 



2 The model and its inputs 

The statistical hadronization model (SHM) [TUB] assumes that all heavy quarks (charm 
and bottom) are produced in primary hard collisions and that their total number stays 
constant until hadronization. The validity of this latter point is investigated in the next 
section. Another important factor is thermal (but not chemical) equilibration in the QGP, 
at least near the critical temperature, T c . The quarkonia are then all produced (non- 
perturbatively) when T c is reached and the system hadronizes. Our picture is different 
from the one of color screening [T7] (see [18] for a more recent account), which assumes 
that formation of quarkonia takes place before a thermalized QGP is established but is 
suppressed in the QGP when the temperature reaches a certain threshold, which is species 
dependent. Recent results from solving Quantum Chromodynamics (QCD) on the lattice 
indicate indeed that J/ip mesons could survive in the QGP up to 1.6T C [15], implying, e.g. 
little melting at SPS energies. Further investigations indicate that T mesons may survive 
up to temperatures of at least 2.3T C [20] . 

In the SHM, we assume that (i) no quarkonia production takes place before the formation 
time tq of the QGP or (ii) that all quarkonia formed before tq are melted in the initial hot 



QGP phase. Under this scenario the possible existence of bound quarkonia states in the 
QGP is not important, since the estimates made in the next section demonstrate that, 
unless there is a resonance-like enhancement of the cross section, formation of quarkonia 
from uncorrelated c and c quarks in the expanding plasma is very unlikely even if such 
states exist as bound states. 

Thermal equilibration of charm quarks in a QGP is currently investigated vigorously 
both experimentally and theoretically. In particular, elliptic flow of charm quarks [22] has 
been demonstrated within a coalescence approach [23"f2l] to be a very sensitive probe 
of thermalization. Large elliptic flow is obtained in hydrodynamic simulations based on 
the diffusion coefficient of charm quarks in a QGP [22]. Elliptic flow of J/ip is smaller in 
recent hydrodynamic calculations [26] than in coalescence models [27] ■ Remarkably, the 
experimental results on the elliptic flow of electrons from charm mesons ("non-photonic" 
electrons) [28.29J do indicate that charm quarks indeed thermalize to a significant degree. 
In addition, a large energy loss of charm quarks in the QGP [30] has been recently found by 
experiments [2Uf 311132] . The experimental results of PHENIX [31 J indicate in fact a larger 
energy loss than theoretically expected [33]. The mechanism for charm thermalization is 
not well understood. A model based on resonant rescattering of heavy quarks in the QGP 
[31] reproduces quite well the experimental results. It is worth noticing that in this model 
the thermalization time for heavy quarks is substantially reduced compared to the case 
of perturbative interactions [31]. In transport models the measured elliptic flow can be 
reproduced only for a charm quark scattering cross section that is much larger than the 
cross sections computed within perturbative QCD [351136] . In any case, we note that the 
SHM can be applied even if charm thermalization is reached only close to T c . 

In the following we briefly outline the calculation steps in our model [TH5] . The heavy 
quark (charm or bottom) balance equation PQ, which has to include canonical suppression 
factors whenever the number of charm or bottom pairs is not much larger than 1, is used 
to determine a fugacity factor g c via: 
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Here N^ r is the number of directly produced cc pairs and /„ are modified Bessel func- 
tions. In the fireball of volume V the total number of open (iV*^ = n™.V) and hidden 
(iV*g = n^V) charm hadrons is computed from their grand-canonical densities and 
n*g, respectively. The densities of different particle species in the grand canonical ensemble 
are calculated following the statistical model [3711381139] . The balance equation (CQ) defines 
the fugacity parameter g c that accounts for deviations of heavy quark multiplicity from 
the value that is expected in complete chemical equilibrium. Since in the SHM the charm 
and bottom quark yields are completey determined by their production in initial (hard) 
collisions the fugacity factor can be very large, especially for bottom production. The 
yield of charmonia of type j is obtained as: Nj = g^.Nj . 

To illustrate some essential features of the model, we show in Fig. [1] the centrality de- 

1 The cross sections for production of charmed quarks are much too small to allow for their 
chemical equilibration in a QGP at reasonable (T <1 GeV) temperatures, as is demonstrated 
in [H]. 




Fig. 1. Centrality dependence of the canonical suppression for charm, (left panel) and of 

the charm quark fugacity, g c (right panel) for SPS, RHIC and LHC energies. 

pendence of the canonical suppression factor, h/Io, and of the charm quark fugacity, g c , 
for SPS, RHIC and LHC energies. The inputs for the calculations are discussed below. 
It is clear that at SPS and RHIC energies, where the number of c c pairs especially for 
more peripheral collisions is not much larger than 1, the centrality dependence of the 
J I ip yield, controlled by the g c parameter, is determined by the canonical suppression. At 
LHC energy this is a rather weak effect, owing to a much larger charm production cross 
section and to a larger volume at freeze-out (hadronization) . We note that an important 
check of the validity of the statistical hadronization model is the if)'/J/i[) yield ratio. The 
model predictions agree with the earlier data at the SPS energy [Tj and we will discuss 
the present status in Section HI The model has the following physical input parameters: i) 
the temperature, T, and baryochemical potential, //&, for statistical model calculations; ii) 
the heavy quark production cross section in nucleon-nucleon interactions; iii) the volume 
of one unit of rapidity, Va v =i at chemical freeze-out. We discuss them below, with special 
emphasis on the updates used here relative to our previous set of predictions [5]. 

Based on our recent analysis of hadron yields within the thermal model [39J , we use here 
(T,fi b ) of (161,22.4) MeV for RHIC and (161,0.8) MeV for LHC. The uncertainty in T 
is 4 MeV and we have shown previously j5] that the outcome of SHM results exhibit 
little sensitivity to the precise value of T. Note that the thermal parameters have been 
determined from fits of data at midrapidity in central collisions. From the measured 
particle abundancies [40], the centrality dependence of T and /x& at RHIC is expected 
to be weak. Based on the measured rapidity dependence of particle ratios at RHIC [43] 
it is expected that will be different away from midrapidity (in particular, //& is 

expected to increase), but this has a rather small influence on the model results in the 
present context. 

In order to have accurate calculations, it is important to include the complete spec- 
trum of open charm states. Based on the latest updates of the Particle Data Group 
[41] . compared to our previous results [5], we have added for the present calculations 10 
new charmed mesons (D^(2400), L'i(2420) ± , £>* (2317), £> s i(2460)) and 12 new charmed 



baryons (£ c (2520), H°(2646), H^ + '°) plus their anti-particles. We now have in the code 
(including antiparticles) 40 charmed meson states and 32 charmed baryon states. As a 
consequence of these updates, the J/0 yield is reduced by 12% at RHIC and by 16% at 
LHC energy as compared to the results obtained with the previous, less complete, spec- 
trum. Even if the spectrum of charmed hadrons may still not be complete at this stage, 
the effect of further missing states on the J/0 yield is suppressed because of their mass. 

To get a feeling for the size of the effect we have investigated the following scenario: 
some of the above new meson resonances are interpreted [42] as chiral partners of the 
known charmed hadrons. Within this approach, their mass is about 300 to 400 MeV 
larger than for the "ordinary" charmed hadrons. It is easy to estimate that, in a thermal 
approach, the abundancy of all possible still missing charm states is about 10% of the 
abundancy of "ordinary" charmed hadrons, which would be equivalent to a reduction of 
the charm production cross section in proportion. As we discuss below, the uncertainty in 
the charm production cross section is much larger and dominates the uncertainty of our 
model predictions. 

The spectrum of open bottom mesons and baryons is poorly known in comparison to 
the charm sector. To partially overcome this difficulty, we have added, besides the known 
states [41] , excited states with masses derived via an analogy with the charm mesons and 
baryons. Despite this, we should clearly state that our predictions on T production, which 
we pursue here for the LHC energy, are an upper bound. 




Fig. 2. Rapidity dependence of the charm production cross section in pp collisions. Left panel: 
for RHIC energy, pQCD calculations [45J and experimental values |50|51| . Right panel: for LHC 
energy, NLO calculations of ref. [M| scaled up by a factor of 1.5. 

The charm production cross section in nucleon-nucleon interactions is, in absence of direct 
measurements, taken from perturbative QCD (pQCD) calculations. For SPS energy we 
use the rapidity density of the charm cross section of ref. |44j . For the RHIC energy, we 
have adopted the recent pQCD calculations [45], which include a careful estimate of the 
theoretical systematic errors. It was shown earlier [44] that, at RHIC energy, shadowing of 
the parton distribution functions (PDF) in AA collisions is a minor effect. For LHC energy 



we have adopted the NLO calculations of ref. [44], where shadowing for AA collisions is 
included. In these calculations the intrinsic kj< broadening of partons ((k^) = l GeV 2 ) is 
supplemented by a nuclear broadening of 0.7 GeV 2 (A=200). Further, we have scaled up 
these pQCD calculations by a factor of 1.5 to bring them in agreement with more recent 
results [46] (see discussion in [H]). Following the study made in [15] we have assumed the 
systematic errors of the cross section to be a factor 2 for SPS and LHC energy. From the 
cross section, we calculate the number of produced cc pairs as: N c5 = a c5 ■ Taa, where Taa 
is the nuclear overlap function, calculated using [4*9] . 

The rapidity distributions of the charm production cross section in pp collisions for RHIC 
and LHC energies are shown in Fig. [2J As seen in the left panel in Fig. [21 at RHIC 
the extracted experimental charm production cross section at midrapidity of PHENIX 
[5T] of 123 ± 47 fib is, within the errors, in agreement with the pQCD calculations of 
da^ CD /dy = 63.7^ j fib [45J. The measured value of 300 ± 98 fib by STAR [50] is 
significantly larger. Note that the experimental values are extracted mostly from single 
electron spectra [50f5"T] , a difficult measurement which is reflected in the size of the errors. 
It was shown [15] that the PHENIX single electron transverse momentum (p t ) spectra [52] 
are compatible with the upper limit of the systematic errors for the pQCD calculations, 
as also seen in Fig. [2] for the pt-integrated value. We note that, at Tevatron, the measured 
cross section of D mesons [53] was underpredicted by pQCD calculations [54J by up to a 
factor of 2 for the lowest measured transverse momentum (~ 5 GeV/c). However, in more 
recent calculations [55] this discrepancy is reduced to a factor of 1.5, so that now data 
and theory are compatible within errors. 




Fig. 3. Rapidity (left panel, for central collisions) and centrality dependence (right panel) of the 
number of bb pairs per unit of rapidity. The lines are for the central value of the bottom cross 
section [H], the shaded areas are bounds for the uncertainty. 

For bottom production cross sections we use the pQCD calculations of ref. [44] . For the 
uncertainty of the cross section we assume a factor of 1.5 up and down. In Fig. [3] we 
show the bottom quark pair production rates as a function of rapidity for central Pb-Pb 
collisions and as a function of centrality for y=0 and y=3. 



In our approach all charmonia are formed by cc recombination at the phase transition. An 
interesting question is the volume within which recombination can take place. A natural 
size for this quantity is the volume corresponding to a slice of one unit of rapidity. We 
have shown earlier [5] that, at RHIC, the dependence of the results on the magnitude of 
the rapidity width is small, as long as intervals of 1-3 units of rapidity are considered. 
This dependence is even smaller for the LHC case. The volume Va^=i is obtained from 
the calculated thermal densities and the experimental values of charged particles rapidity 
densities, dN^/dy, at midrapidity. For central collisions (number of participating nucleons 
N part =350) we obtain at SPS Va„=i=1200 fm 3 , while at RHIC, dN ch /dy=701 leads 
to VAy=i=2400 fm 3 [3H]. At LHC, based on a recently-proposed parametrization [SB] 
of the energy dependence of dN c h/dy, the expected value is dN c h/dy=1816, leading to 
Va^=i=6200 fm 3 at midrapidity. This presently used value of dN c h/dy at LHC is slightly 
lower than our previous phenomenological extrapolation [5] . We have shown earlier [5] that 
the uncertainty in the volume is affecting the model predictions rather marginally. Our 
previous results [5] were obtained assuming a linear scaling of dN ch /dy (and consequently 
Vaj/=i) with Npart. For the present calculations the scaling is done for the 'core' region 
only (see below). 



Fig. 4. The nuclear density (Woods-Saxon) distribution for the Au nucleus. The core corresponds 
to the radius Ra + X c , beyond this radius the nucleons are in the corona region. 

An important effect to be included when studying the centrality dependence in nucleus- 
nucleus collisions is the so-called "corona effect" . Nucleons from the surface of the colliding 
nuclei do not participate in the formation of the hot fireball (core) were QGP is assumed to 
be produced. This fact has recently been emphasized also in ref. [57]. Since the SHM applies 
only to the QGP zone, and since charmonium production in nucleon-nucleon collisions 
is, in general, very different from that predicted in the SHM, it is relevant to distinguish 
between core and corona in our approach. To quantify core and corona, we use calculations 
of the nuclear overlap \4^ , employing a Woods-Saxon nuclear density distribution, plotted 
for the Au nucleus in Fig. ^ We define the core region as corresponding to the half-density 
nuclear (charge) radius (Ra=6.37 fm for the Au nucleus) extended by a thickness X c (see 
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Fig. and derive in this manner N part and N co u for the core and corona regions. The 
fraction of participating nucleons contained in the corona region is shown as a function of 
centrality in the left panel of Fig. [3]/br three values of X c . Effectively this three-dimensional 
approach contains a sharp transition between core and corona. We treat the core as QGp 
using the SHM and the corona as superposition of nucleon-nucleon collisions. 




Fig. 5. Left panel: the fraction of nucleons in corona as a function of centrality. Right panel: the 
centrality dependence of charged particle yields at midrapidity. The data [59] are compared to 
estimates based on the corona contribution (the curves are arbitrarily normalized). 

It was recently demonstrated [57] that the centrality dependence of particle production at 
RHIC and SPS can be well explained assuming an interplay of bulk production in the core 
and of elementary collisions in the corona part. Within this picture, we compare in the 
right panel of Fig.[5]our expectation for the centrality dependence of charged particle mul- 
tiplicity with RHIC data from PHOBOS [58], PHENIX [59] and BRAHMS [60]. For the 
elementary collisions we use diV c / l /d77=2.5 (at midrapidity, non single-diffractive), mea- 
sured in pp collisions at -^5=200 GeV [BT] . We normalize our calculations to the measured 
value [59] for A^ art =300 to extract for the core dN c h/drj/ (N part /2)= 4.23, 3.92 and 3.81 for 
X c = 0, 0.7 and 1.2 fm, respectively. We conclude that this treatment of the corona effect 
reproduces the trend in the data only for the case X c =0 which we consider somewhat 
extreme since, effectively, pA collisions do not contribute to core particle production (see 
above). We note also that, using BRAHMS data, one may come to a somewhat different 
conclusion. We consider for the corona thickness A c =1.2 fm as a realistic value, since at 
that point the density has dropped to about 10 % of the central value and adopt it as a 
baseline choice for our following model calculations. In this case our calculations account 
for the measured centrality dependence of bulk hadron production only partially. On the 
other hand, the interpretation of such data within saturation ("color glass condensate") 
models [62 56] is quite successful. However, it is clear from the present exercise that the 
corona contribution needs to be considered in order to have a realistic modelling of the 
collision [57]. To let the reader judge the influence of the corona effect we provide calcu- 
lations for the centrality dependence of J/ip production with and without considering the 
corona. 




Fig. 6. Centrality dependence of rapidity den- 
sity of charm yields at RHIC and LHC ener- 
gies for the total overlap (open symbols) and 
for core only (full symbols). 



Fig. 7. Energy dependence of rapidity density 
of J/ip cross section in pp(p) collisions (see 
text for details). The LHC energy for Pb+Pb 
is marked by the vertical line. 



Based on the central values of the pQCD charm production cross sections |4"5]I44"] . we show 
in Fig. [6] the centrality dependence of the number of cc pairs for central rapidity for the 
RHIC and LHC energies. In the following we perform calculations within the statistical 
hadronization model for the core, and add the contribution of the corona region for which 
we use the J /if} cross section in elementary collisions. The energy dependence of the J /if) 
cross section (da/dy) in pp collisions is shown in Fig. [7J Recent experimental data at RHIC 
from PHENIX [63] and at the Tevatron from CDF [64] are compared to calculations by 
Gavai et al. [65] for two PDFs. In our model calculations we employ for RHIC the measured 
[63] J /if) cross section of 0.724 fib, which is about 1% of the pQCD charm production 
cross section [43]. Assuming the same 1% of the calculated pQCD charm cross, for LHC 
we derive a J ' /if) cross section of 6.39 fib, which we use in the following calculations. This 
cross section is in line with the calculations [65] which reproduce the Tevatron data [64] 
(see Fig. [7j). For SPS energy we employ a J/ 'if) cross section at midrapidity in pp collisions 
of 50 nb, obtained from an interpolation based on a recent compilation of data [66]. These 
values are summarized in Table [TJ 



Table 1 

Summary of the values used in our calculations for charm. The values for da^/dy are from 
pQCD calculations, the values for da p Jj^ are from measurements, either interpolated (for the 
SPS energy, cf. ref. [66]), or directly measured at RHIC [63J or assumed as 1% of da^/dy at 
LHC energy. 
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d^£/dy (jib) 


daf^/dy ( M b) 


17.3 




0.050±0.030 


200 


63.7±H;| 


0.774±0.124 


5500 


639lSg 


6.4±3.2 



The yield of J/ip mesons, e.g., is then calculated as follows: 
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and we use for of% el values of 30, 42, and 60 mb at SPS, RHIC, and LHC energies, 
respectively. 

The T(1S) cross section in pp collisions at the LHC energy is assumed to be 36.8 nb, 
which is 1.7-10~ 3 of the bb cross section at v /s/v r /v=5.5 TeV. This is the fraction derived 
from the Tevatron data at 1.8 TeV [611671168] . 
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Fig. 8. Volume dependence of the rapidity density of the J/tp yield calculated within the SHM, 
normalized to the cc yield. The measured value for the RHIC energy }63j51j is indicated by the 
shaded band, while the open box depicts the interpolated value for the SPS energy. 

Before proceeding to discuss our calculations for nucleus-nucleus collisions, we investigate 
the model predictions for the proton-proton case. In Fig.[8]we show the volume dependence 
of the predicted J/ip yield normalized to the cc yield. It is seen that, for a realistic 
volume expected for pp collisions, of about 30-50 fm 3 , the model strongly overestimates the 



measured [63"f5"T] or interpolated [66] values for the RHIC and SPS energies, respectively. 
This is expected, as charm quarks are not likely to thermalize for small volumes. As 
a consequence, we consider, somewhat schematically, the volume for which the model 
calculations reach the measured values in pp collisions to be the minimal core volume 
required for the validity of our model. Based on the results shown in Fig. [HI we adopt 
the value Vggp=400 fm 3 . Below this value, we consider that the J/ip production takes 
place exclusively via hard processes in elementary collisions. In nature, the transition from 
the QGP region to the pp region is likely to be continuous, but we have no quantitative 
means to model this transition. Also, from Fig. [S] it is clear that Vqqp carries a rather 
large uncertainty. In addition, we do not consider the experimental fluctuations in the 
centrality selection which would further smoothen the transition. 



3 Annihilation of charm quarks in the QGP 

In this section we investigate the assumption, central to our model, that the number of c 
and c quarks which are produced in initial hard collisions stays constant during the evo- 
lution of the QGP until the critical temperature T c is reached and the system hadronizes. 
For such an estimate we assume that the light quarks and gluons reach equilibrium to 
form a QGP relatively quickly after the first hard collisions, on a time scale of the order 
of 1 or at most a few fm/c. We first note that further production of cc pairs in the QGP 
can be completely neglected at temperatures T < 1 GeV |21j . as the charm quarks are far 
from chemical equilibrium. On the other hand, charm quark annihilation may take place 
via reactions of the type 

c + c-^g + g (5) 



or via more complex reactions such as those forming 3 or more gluons, a quark- ant iquark 
pair plus gluons, etc. Within the framework of perturbative QCD, the cross sections for 
Eqs. [51 [6] can be evaluated. In the following, we perform such an evaluation for Eq. [5] 
and use it to determine the annihilation rate of charm quark pairs in the QGP via this 
channel. The 2-gluon annihilation is the dominant channel, its cross section significantly 
exceeding all other listed cross sections [61] . 

Channels with 3 or more gluons are not expected to play a major role in our estimate. An 
estimate of the relative size of annihilation into 3 gluons vs. 2 gluons can be obtained by 
inspecting the width of the two lowest-lying charmonia r] c and J/ip. Because of C-parity 
conservation the r] c 's width is dominated by 2-gluon decay while J/^> decays mainly into 
3 gluons. In fact, their measured widths are 25 and less than 0.1 MeV, respectively. To 
leading order in a s and in non-relativistic approximation, the ratio of these widths is: 



or 



c + c 



q + q, 



(6) 



r(J/V) _ 4 



(tt 2 - 9)a s . 



(7) 



Evaluation of this expression for a s (m(J/ip) 2 ) ~ 0.3 [70] overestimates the measured ratio 
by approximately an order of magnitude. Annihilation into 3 gluons is obviously strongly 
suppressed relative to the 2-gluon process. For more information on this issue see [71] . 
Moreover, multi-gluon annihilation should be further suppressed by phase space if gluons 
acquire, as expected, a thermal mass in the QGP. We conclude that annihilation into 2 
gluons is the dominant channel of the inelastic cc cross section. 



To evaluate the total annihilation rate we start from the rate equation 



dr c5 
dr 



(8) 



where (a C c^ gg v r ) is the thermal average of the annihilation cross section times the relative 
velocity v r in the QGP, and n c = n e is the charm quark density. The quantity is 
the annihilation rate per volume or the rate of change of the charm quark density. The 
cross section for Eq. [5] is taken to first order in the strong coupling constant a s from the 
work of |72j, where the inverse cross section for g + g — > c + c is computed, applying 
detailed balance. To get an estimate of the upper limit for this cross section in the QGP 
we evaluate the equation of [72] using a s = 1 and for a charm quark mass of m c = 1.5 
GeV. 

The charm quark density n c depends on the evolution time of the QGP, i.e. 

dNjdyjr) dN c /dy(r ) 
Hc V(Ay = l,r)-V(Ay=l,ry [J) 



where tq is the initial time of QGP formation and dN c /dy is the charm quark rapidity 
density. 

The total annihilation yield of charm quarks in the QGP is then given by 

/ c dr - 
-^V(Ay = l,r)dr. (10) 

TO 



To proceed further we need to model the evolution of the QGP. In the spirit of our 
upper limit estimate we assume here a 1-dimensional Bjorken-type expansion of the QGP 
yielding, using entropy conservation, a relation between T and r: 

^(32 + 21^)^ = 3.8^, (11) 
45 v A\ 



where Nf is the (effective) number of massless flavors, and dN/dy is the total particle 
rapidity density. The transverse system size A± is about 150 fm 2 at T c for central Au-Au 
or Pb-Pb collisions. Using this scenario we get for the temporal evolution of the volume 
then 



V(Ay = l,r) = A ± r. 



(12) 



The total annihilation yield of charm quarks in the QGP is, using Eqs. I9fl0lll2l given by 



'dN r , A 2 1 fdr 



h^W -J-J-^Vr). (13) 



T"0 



In the following we consider 2 scenarios (assuming Ny = 2.2 and r = 1 fm for both): i) at 
RHIC energy, this leads to r c = 2.7 fm and to an initial temperature T(r ) = 225 MeV for a 
charged particle rapidity density of 660 and ii) the LHC energy scenario results in r c = 8.3 
fm and T(ro) = 325 MeV for a charged particle rapidity density of 2000. Following [73] we 
evaluate the temperature dependence of the thermal average in Boltzmann approximation. 
The resulting temperature dependence of the thermal average is presented in Fig. EH 
Since the charm quark annihilation cross section decreases with increasing scale (\/s), the 
resulting thermal average also drops with increasing temperature. 

To get numerical results on the annihilation rate we determine the integral in Eq. [13] for 
both scenarios numerically. Note that the lifetime of the QGP enters (to 1st order) only 
logarithmically, implying that the details of the expansion are not very important for our 
estimate. In particular, using much smaller values for r and correspondingly higher initial 
temperatures changes the results only marginally. We furthermore add that the volume 
used here for an upper limit on the annihilation yield is computed when the system 
reaches T c . Taking into account the increase of volume until chemical freeze-out, which is 
used in the next section, would further decrease the annihilation yield and, hence, is not 
considered here. 
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Fig. 9. Temperature dependence of the ther- Fig. 10. Annihilation rate as a function of ini- 
mal average as defined in Eq. [8] (see text). tial charm rapidity density. 

The final results of the estimate for the fraction of charm quarks lost via annihilation, 
5 = N^ m j dN c / dy(r ) are given in Fig. [10] for the RHIC and LHC scenarios. Since the 
annihilation yields increase quadratically with the charm quark rapidity density we start, 
in the spirit of an upper limit estimate, the calculation at the values predicted by pQCD 



calculations for each scenario, i.e. -^-ijo) = — ^ — = 1-4 and 13.3, and then increase 
these numbers by up to a factor 5. 

In general, the fraction of annihilating cc pairs in the QGP is negligible for both sce- 
narios. Even for an unrealistically high value of ^-(to) ~ 60 in the LHC scenario, the 
total fraction of annihilated pairs is less than 1%, and decreases linearly for lower charm 
quark densities. For realistic scenarios our estimate implies that charm quark annihila- 
tion in the plasma can be safely neglected. Along the same line, production of charmonia 
via uncorrelated charm quark annihilation in the QGP is expected to fall significantly 
below the above computed annihilation yield into gluons, lending strong support to our 
interpretation that all quarkonia are produced late, when the system reaches the critical 
temperature and hadronizes. 

This result has a number of significant consequences. We first note that the small anni- 
hilation rate in the plasma also implies that charmonia are not likely to be formed in the 
plasma. Therefore, the possible existence of (quasi-)bound J /if) mesons in the QGP [741T9] 
would not influence our predictions for charmonia yields as these states would never be 
populated during the evolution of the QGP. We neglect in this context J /if) production 
in hard collisions before QGP formation, i.e. before r , because of the the large formation 
time [75] . 

These considerations also highlight the differences between our approach, where all char- 
monia are formed non-perturbatively at T c , during hadronization of the QGP, and the 
recombination model of [7118] , where the cross section for the production of J /if) mesons 
from c and c quarks is obtained by first computing, using the operator product expansion 
technique [7B], the dissociation of J /if) mesons due to collisions with gluons and then in- 
verting this cross section using detailed balance. We note that this procedure yields a cross 
section for J / if) production of several mb [9] , significantly exceeding that for annihilation 
of the cc pair into 2 gluons, see Eq. [51 



4 Centrality and rapidity dependence of quarkonia yields 

In our earlier studies [Tf5] we have shown that the J /if) data at SPS energy can be described 
within the statistical approach, but only when assuming that the charm production cross 
section is enhanced by about a factor of 3 beyond the perturbative QCD (pQCD) predic- 
tions. Those results were obtained assuming that hadronization of charm takes place in 
the whole volume of the system at chemical freezeout. We consider this assumption to be 
rather extreme and adopt here a picture of hadronization within one unit of rapidity at 
midrapidity, as done for the analysis of the yields of hadrons without charm within the 
thermal model [37|39] . The comparison of our model calculations with the data measured 
at SPS by the NA50 experiment is shown in Fig. [TTJ Two sets of data are included in 
Fig. dU the NA50 data of 1998, as analyzed by Gosset et al. [77J and with a further 
normalization (see pQ) and the NA50 data of 2004 [78J, representing the yield of J /if) nor- 
malized to the Drell-Yan yield, which we arbitrarily normalize for the sake of comparison, 
assuming that the Drell-yan yield scales with N co u. It is seen that the two data sets show 
a very good consistency in the centrality dependence. Note that a factor of 0.5 was used 



to convert the full phase space experimental data [771fT] . to midrapidity yields. 



The experimental data are well described by invoking a moderate enhancement of the 
charm production cross section of a factor of 2 compared to pQCD calculations [2]. Very 
clearly, adopting for the present calculations the hadronization within one unit of rapidity 
leads to an agreement with the data for a smaller enhancement of the charm cross section 
than required in our earlier study [5]. New measurements by the NA60 experiment [81] 
indicate that the enhancement in the dimuon yield below the J / ip mass, earlier observed by 
the NA50 jE2] , is of prompt origin and, as such, cannot be interpreted as an enhancement 
of the charm production cross section, although an experimental result on the charm cross 
section is currently not available. We also note that a factor of 2 enhancement is within 
uncertainties of the pQCD calculations. 




Fig. 11. Centrality dependence of J ftp yield per number of collisions (left panel) and of the ratio 
of yields of and J/ip (right panel) at the SPS energy. The shaded band in the left panel 
denotes the range in the charm production cross section (a factor of 2 up and down) and the 
error assumed for the J ftp cross section in pp collisions. The experimental data on J/ift yield are 
from ref. [77] (see [I]) and ref. [78] (2004 NA50 data) and on the ip'/J/ip ratio from ref. |79IS()j . 

A comparison between calculations and data [79|l80] is shown in the right panel of Fig. [TTJ 
for the centrality dependence of the yield of ip' relative to J/tp. Using the present core and 
corona model leads to a 25% increase of the ratio compared to the case when all nucleons 
are considered to be part of the QGP region (no corona). The model predictions are shown 
only down to the N part values corresponding to Vqgp- For more peripheral collisions the 
data suggest a continuous transiton towards the value measured in elementary collisions. 

For RHIC energy, the centrality dependence of the J ftp rapidity density at midrapidity 
is shown in Fig. [TJl considering three cases for the charm production cross section: i) as 
calculated in pQCD 05], and as measured by ii) PHENIX [BTJ and iii) STAR [50J. The 
model calculations agree with the recent PHENIX data [H] very well for the central value 
of the pQCD charm production cross section, in line with our earlier comparison [5]. The 
yield is scaled by the number of binary collisions, N co u. In this scaled representation the 
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Fig. 12. Centrality dependence of J/ip 
yield at midrapidity per number of colli- 
sions at RHIC. The experimental data are 
compared to three cases for the model cal- 
culations, considering the charm produc- 
tion cross section: i) as calculated in pQCD 
05], and as measured by ii) PHENIX [5]] 
and iii) STAR [50] experiments. Note the 
different scales on the vertical axis. The 
lines are our calculations for the central 
value of the charm production cross sec- 
tion, the band corresponds to its system- 
atic uncertainties. The dashed lines are 
model predictions for the rapidity y=1.7. 
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model predicts a slight decrease of the J/ip yield as a function of N part , a trend which is 
compatible with the data. The model overestimates the data by about a factor of 2 for 
the central value of the charm production cross section measured by PHENIX |51j, but 
is compatible with the data for the lower limit of this measured cross section. If we use 
the values of the charm cross section measured by STAR [50] the model overestimates the 
data by about a factor of 10 and is clearly incompatible with the data also concerning the 
trend function of centrality. 

Note that the steps in the calculations seen around N part =70 are the outcome of the 
sharp transition from the core region to the pp region below Vqqp=400 fm 3 . As we have 
mentioned earlier, this transition is probably continuous, but we have not tried at this 
stage to model it. 

The rapidity dependence of the J / ip yield is shown in Fig. [13] for two centrality bins for 
Au-Au collisions. The PHENIX data [H] are well described by the model calculations 




Fig. 13. Rapidity dependence of the J/ip yield at RHIC for two centrality bins. The SHM 
calculations, performed for the nominal pQCD charm production cross section (continuous line 
with the band denoting the systematic errors of the cross section) are compared to recent 
PHENIX data [H]. 

for the central value of the pQCD charm production cross section. In particular, the 
narrow rapidity distributions predicted within the kinetic recombination model [9] are not 
observed here because of the rather broad rapidity distribution of the charm production 
cross section [16]. Based on the PHOBOS measurement of the pseudorapidity density of 
charged particles [83], we have assumed a constant volume Vaj,=i as a function of rapidity. 
As a consequence, the rapidity dependence is determined in our model solely by the input 
charm cross section. 




Fig. 14. Predictions for the J/ip rapidity density at LHC. Left panel: centrality dependence 
(normalized to N co n) at midrapidity, right panel: rapidity dependence for central collisions 
(N part =350). 



The model predictions for the centrality and rapidity dependence of J/ip production 
at LHC are shown in Fig. [TH For most of the range of the charm production cross 
section considered, a growth with centrality faster than that of the number of binary 
collisions is seen is seen in the J / ip yield, allowing to make a distinct case for the statistical 
hadronization scenario at LHC. This increase is stronger for larger charm cross section 
values. The rapidity distribution is broader than at RHIC energy. As a consequence, the 
measured rapidity densities of J/ip in the ALICE central barrel (\r]\ < 0.9) and in the 
muon spectrometer (rj = 2.5 — 4.0) are expected to be comparable. 




Fig. 15. Yield of T at LHC. Left panel: centrality dependence for y=0 (continuous line with error 
band for the uncertainty of bottom cross section) and y=3 (dashed line). Right panel: rapidity 
dependence for central collisions (N part =350). 

The model predictions for the centrality and rapidity dependence of T production at 
LHC are shown in Fig. [151 A decrease as a function of centrality is seen (note again the 
normalization to N co u), similar to the one observed for J/ip at RHIC. 

Fig. US] shows a comparison of the calculated J/ip yield per number of cc pairs for RHIC 
and LHC energies. This figure illustrates the succinct change in trend expected for mea- 
surements at LHC energy, due to the strong increase with energy of the total charm 
production cross section. Note, in particular the decrease with centrality of the J / ip yield 
at RHIC, originating from the canonical suppression of the open charm hadrons. The dot 
corresponds to the PHENIX measurement of J / ip production in pp collisions [63J normal- 
ized to the cc cross section, also measured by PHENIX [51J . This value is smaller than 
the value of 1% derived using the pQCD value [15] for the charm cross section and much 
smaller than the value of 2.5% derived earlier by Gavai et al. |65j . 

In the left panel of Fig. [17] we present the fraction of J / ip and of T yields from the core 
(QGP) region relative to the total (core and corona) overlap. For the J/ip yields, this 
contribution of corona to the total J/ip yield is moderate for all the three energies. Note 
as well the different centrality dependences for the three energy regimes. The T yield is 
predominantly originating from the core part. The right panel of Fig. [T7I shows the factor 
by which the J/ip yield increases when considering the corona contribution, compared 
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Fig. 16. Centrality dependence of J/ip yield normalized to cc yield for the SPS, RHIC and LHC 
energies. The thick line is the value used in the model for pp collisions. The dot corresponds to 
the PHENIX measurements in pp collisions at RHIC energy. 




Fig. 17. Centrality dependence of corona contribution. Left panel: the fraction of yields of J/tp 
and T from the core (QGP) region relative to the total overlap. Right panel: the increase of the 
J/ip yield considering corona, compared to the case when all the overlap volume is QGP. 

to the case when we assume that the entire overlap volume is QGP. The curves reflect 
the differences between the yields calculated in our model and the yields in pp collisions 
(see Fig. [TBI) . Since we consider a rather modest corona fraction in our calculations, the 



increase of the J/ip yield due to the corona is moderate. 



5 Transverse momentum dependence 



We turn now to the analysis of the transverse momentum (p t ) spectrum of J/ip mesons. 
Within our model, the p t distribution of J/ip, as well as of any hadron carrying charm 
(or beauty), is determined by the temperature and transverse expansion velocity, (3, at 
chemical freeze-out. We employ a formula for collective expansion [S3], but consider, to 
keep the input minimal, one average velocity instead of a velocity profile: 



_dN_ 
Pt ■ dpt 



m t ■ I 



p t smhy t 
T 



PtCOshy t 
T 



(14) 



where m t = ym^ +p%, m is the rest mass of the particle, and y t = tanh _1 (/5), with f3 
the collective (transverse) expansion velocity in units of the speed of light. The 1$ and K\ 
are the modified Bessel functions. 
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Fig. 18. Centrality dependence of the average transverse momentum of J/tp mesons in Pb-Pb 
collisions at SPS energy. The data measured by the NA50 experiment [86J are compared to 
calculations using Eq. [HI 

Our picture is similar to that tested with data at the SPS by Gorenstein et al. [85] (first 
proposed by Grandchamp and Rapp [3]). Spectra of both J/ip and ip' mesons and of Q 
hyperons were shown [85] to be described by a (kinetic) temperature equal to the chemical 
freeze-out temperature, while the expansion velocity was found to be substantially lower 
than for the other hadrons. This indeed suggests a kinetic freeze-out of J/if> coincident 
with the chemical freeze-out (at the phase boundary between QGP and hadronic matter 
|39|), with (3 characterizing the transverse expansion in the QGP phase. 



To check the validity of this interpretation, we calculate with Eq. [TH the transverse mo- 
mentum spectrum of J/if> and compare it to the available data. For SPS energy we show in 
Fig. [TS] the centrality dependence of the average transverse momentum of J/ ip measured 
by the NA50 experiment [SB]. Lacking a measured J/ip spectrum in pp collisions, we have 
not included the corona contribution in our calculations. The data are well described by 
the calculations using T=160 MeV and /?=0.22, in agreement with the earlier results [85] . 
Note that the data approach a constant value for central collisions (N part >150), where 
the present model is applicable. We would like to emphasize that the (p t ) data cannot 
constrain T and (3 in a satisfactory way. When fitting the data, local minima are present 
in the \ 2 distribution. Several sets of T and f3 values describe the data equally well, as 
known from fits of distributions of other hadron species [ST]. This problem would be al- 
leviated if high-precision spectra as a function of centrality become available. We note 
that the centrality dependence of the average transverse momentum measured at SPS 
[SB] , Fig. [TSJ is well explained over the whole centrality range assuming a purely hadronic 
scenario [9"U|91ll92j (see discussion in ref. [TB]). 
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Fig. 19. Transverse momentum distribution of J/ip in Au-Au collisions at RHIC energy |14j . 
The data (symbols) are compared with calculations using Eq. [TH for two cases of the kinetic 
freezeout parameters (lines). The spectrum measured in pp collisions [88] with the fit as used 
for the corona part in our calculations is also shown. 

For RHIC energy we calculate the J /if) spectrum taking into account both the QGP con- 
tribution (according to Eq. UM and the corona part. For the corona, we use the spectrum 
measured in pp collisions [88] . which we fit in order to be able to interpolate or extrap- 
olate. Fig. [19] we show the comparison of the J /if) spectra measured by the PHENIX 
experiment [H] to our calculations. Two sets of parameters are investigated for the QGP 
part: i) T=160 MeV, /?=0.3 and ii) T=360 MeV, /?=0.1. The data are compatible with 



the picture of J/ip production via hadronization at chemical freeze-out, expected to be 
described by scenario i). In detail, the data appear to exhibit a high-pt component which 
is not accounted for by the corona component used in our model and which seems better 
described by a larger temperature in the QGP region. Interestingly, the expansion velocity 
extracted from fits of the spectra of Q and particles, which are also expected to interact 
rather weakly in the hadronic stage, yields for T=160 MeV (3 ~0.4-0.5, significantly lower 
than (3 values determined for other hadrons [87]. We note that spectra of open charm 
mesons, as inferred from the measured nuclear modification factor of non-photonic elec- 
trons [2H] are also compatible with a thermal distribution at chemical freeze-out with a 
flow velocity (3 ~0.3 as used here for the J/ijj spectra. 




Fig. 20. Predictions for momentum spectrum of J/ip (left panel) and T (right panel) for different 
values of the expansion velocity, /?, for central Pb+Pb collisions (N part =350). Also included are 
the measured spectra in pp collisions at Tevatron |64J. 

In Fig. [20] we present predictions for the transverse momentum spectrum of J / ip and T 
in central collisions at LHC. We show the spectra for a range of expansion velocities, 
as expected for LHC energy. For pp collisions in the corona we have assumed the same 
spectral shape at LHC as measured in pp collisions at Tevatron [64], which are also 
included for comparison in Fig. [20] The corona fractions corresponding to central Pb+Pb 
collisions (N part =350) are 0.13 for J/ip and 0.09 for T. For both J/ip and T, the spectral 
shapes expected within the statistical hadronization model are very different compared 
those expected from the superposition of pp collisions, making the data at LHC a stringent 
test of our model. 



6 Conclusions 



We have presented new results on the statistical hadronization of heavy quarks at the 
SPS, RHIC and LHC energies. We have focused on the results for J/ip yields, but have 



also extended the model to predict T yields in Pb+Pb collisions at the LHC energy. 
Included in our model is now a separation of the collision geometry into a core (QGP) 
and a corona (pp collisions) part. Our estimate of the annihilation rates of charm quark 
in a hot plasma demonstrates that this is a small effect. An important ingredient in our 
model is the charm production cross section, which is presently quite poorly measured. We 
have considered the effect of its uncertainty on the model predictions. For SPS energy, our 
present calculations are performed assuming that the hadronization takes place within one 
unit of rapidity at midrapidity, as used also for RHIC and LHC energies. Good agreement 
with data is found with a charm cross section which is moderately (a factor of 2) larger, 
but consistent within errors with pQCD calculations. For RHIC energy we have studied 
the centrality and rapidity dependence of J/ip production. Our model is in good agreement 
with the available data, lending strong support for the picture of statistical hadronization. 
We have discussed the transverse momentum distributions of J/if> mesons expected from 
the model. The present data, in particular those at RHIC energy, are compatible with 
the picture of the J/if> (kinematic) freeze-out coincident with chemical freeze-out, but 
more precise data are needed to quantitatively explore this idea. Our predictions for LHC 
energy concerning centrality, rapidity and transverse momentum dependence of the J/ip 
yield will be tested with data within the next few years. 
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